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ABSTRACT 


Cebeci's  viscous/inviscid  interaction  program  was  applied 
to  the  analysis  of  steady,  two  dimensional,  incompressible 
flow  past  four  airfoils,  the  NACA  663-018,  0010  (Modified), 
4412  and  the  Wortmann  FX  63-137.  Detailed  comparisons  with 
the  available  experimental  results  show  that  the  essential 
features  are  correctly  modelled,  but  that  significant 
discrepancies  are  found  in  regions  of  flow  separations. 
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I .  INTRODUCTION 


Understanding  of  the  characteristics  of  the  airflow  over 
an  airfoil  is  of  paramount  importance  to  the  airfoil 
designer.  Two  methods  are  currently  available  which  give 
accurate  results.  The  first  is  the  use  of  wind  tunnel  tests. 
The  drawbacks  to  this  method  are  cost  and  time  consumption. 
The  second  is  the  processing  of  the  Navier-Stokes  equations. 
This  method's  drawbacks  are  the  requirements  and  expense  of 
using  supercomputers  due  to  the  extensive  calculations  and 
storage  requirements.  There  is  still  a  need  to  come  up  with 
an  inexpensive,  fast  and  accurate  engineering  tool  to  compute 
airfoil  flows. 

Several  methods  have  been  derived  to  accomplish  this  end. 
But  the  most  promising  is  the  Viscous-Inviscid  Interaction 
method.  The  outer  flow  is  computed  using  inviscid  flow 
equations,  and  the  inner  flow  is  computed  using  Prandtl's 
boundary  layer  equations.  The  key  to  this  method  is  the 
extent  of  interaction  between  the  inner  and  outer  flows. 

The  purpose  of  this  thesis  is  to  evaluate  the  capability 
of  the  viscous-inviscid  interactive  aircode  developed  by 
Tuncer  Cebeci  and  associates  at  the  Douglas  Aircraft  Company 
[Ref.  1].  This  computer  program  was  applied  to  four  airfoils 
with  various  angles  of  attack  and  Reynolds  numbers.  The 
computer  results  were  then  compared  to  previously  reported 
experimental  results. 
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The  conservation  of  mass  and  momentum  are  summarized  in 
Chapter  2,  inviscid  flow  calculations  are  discussed  in 
Chapter  3,  and  viscous  flow  equations  are  described  in 
Chapter  4.  Viscous  calculations  are  presented  in  Chapter  5, 
and  the  specific  interaction  methods  are  shown  in  Chapter  6. 
Finally,  in  Chapter  7  computer  and  experimental  results  are 
compared  for  the  NACA  663-018,  0010  (Modified)  and  4412 
airfoils  as  well  as  the  Wortmann  FX  63-137  airfoil. 
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II .  FUNDAMENTAL  EQUATIONS 


The  conservation  of  mass  and  conservation  of  momentum 
provide  the  foundation  for  incompressible  flow  analysis. 
With  these  fundamental  concepts  along  with  appropriate 
assumptions  and  approximations,  working  relations  for  two- 
dimensional,  incompressible  flow  are  obtained. 

A.  CONSERVATION  OF  MASS  (CONTINUITY) 

The  conservation  of  mass  principle  states  that  mass 
cannot  be  created  nor  destroyed.  Equating  this  statement  to 
a  fixed  control  volume  the  net  mass  flow  rate  into  and  out  of 
the  control  volume  equals  the  time  rate  of  change  of  mass 
within  the  control  volume  [Ref.  2:p.  A-l]. 

Given  a  control  volume,  the  mass  flow  rate  through  one  of 
its  surfaces  is  equal  to  the  product  of  the  fluid  density, 
the  fluid  velocity  normal  to  the  surface  and  the  area  of  the 
surface  [Ref.  3:p.  29]. 

d  mass 

-  =  V  •  n  s  (2.1) 

dt 

In  2-D  flow  the  x-component  of  the  mass  flow  rate  at  the 
center  of  the  positive  x-face,  position  dx/2  and  side  length 
dy,  is  represented  by  Taylor  series  expansion  [Ref.  2:p.  A-2] 
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d  dx  02  dx  l 

pu  +  — (pu) —  +  - (pu)(  —  )2  —  +. 

d*  2  Ox2  2  2! 


dy. 


(2.2) 


As  dx  approaches  zero  all  higher  order  terms  disappear 
leaving 


pu 


3  dx 

+  — (pu)  — 

dx  2 


dy. 


(2.3) 


Similarly,  the  x-component  of  the  mass  flow  rate  at  the 
center  of  the  negative  x-face,  position  -dx/2  and  side  length 
dy,  is  represented  by 


pu 


d  dx 
—  ( pu)  — 
dx  2 


dy. 


(2.4) 


As  illustrated  in  Figure  2.1  the  net  mass  flow  through  the 
four  sides  of  the  2-D  control  surface  is 


3 ,  ,dx 

d  dx 

pu  -  — ( pu) - 

dy  - 

pu  +  — (pu) - 

dx  2 

dx  2 

d  dy 

d  dy 

pv  -  — (pv) - 

dx  - 

pv  +  — (  pv) - 

dy  2 

dy  2 

4 


Figure  2.1  Mass  Flow  Through  2-D  Control  Surface 
[Ref.  4 : p .  12] 
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which  is  equal  to  the  rate  of  change  of  the  mass  within  the 
control  volume 


—  dxdy  (2.6) 

3t 

Combining  (2.5)  with  (2.6)  and  simplifying  yields 

-  — (pu)dxdy  -  — (pu)dxdy  =  —  dxdy  (2.7) 

c)x  e)y  c)x 

Dividing  by  dxdy  and  rearranging  yields 

+  JL(pu)  +  <L(pu)  =  0  (2.8) 

c)t  Bx  By 

For  steady,  incompressible  flows  the  continuity  equation 
becomes 


Bu  Bv 

—  +  —  =0  (2.9) 

c)x  By 

or  in  vector  form  the  continuity  equation  [Ref.  3:p.  30]  is 

7  •  V  =  0  (2.10) 

B.  CONSERVATION  OF  MOMENTUM  ( NAVIER-STOKES ) 

The  conservation  of  momentum,  Newton's  second  law  of 
motion,  states  that  the  rate  of  change  of  the  linear  momentum 
is  equal  to  the  sum  of  the  forces  applied  [Ref.  2:p.  B - 1 ]  . 
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(2.11) 


d 

2F  =  — (mV) 
dx 

As  illustrated  in  Figure  2.2  the  two  significant  forces  which 
act  on  an  element  of  fluid  are  surface  forces  which  act  on 
the  surface  only,  pressure  and  shear,  and  body  forces  which 
affect  the  mass  of  the  element,  such  as  gravity.  Assuming 
moment  equilibrium  in  an  element,  xxy  =  xyx,  the  2-D  first 
order  Taylor  series  expansion  for  normal  and  shear  surface 
forces  in  the  x-direction  is 


dx 


dx 


+  “ (  X  x  X  ) ”  X  xx  +  (Xxx) 

dx  2  dx  2 


dy 


3  ,  ,d* 

+  - (Xxy) -  ~  Xxy 

dy  2 


d  dy 

+  ~ (Xxy)  — 

dy  2 


dx 


pj  c) 

—  (Xxx)dxdy  +  — -(Xxy)dxdy. 

dx  dy 


(2.12) 


The  body  forces  per  unit  mass  are  represented  by 


=  Xi  +  Yj  *  Zk 


(2.13) 


such  that  the  x-component  of  the  body  force  on  an  element  is 


f  x  <  body  >  =  pdxdy •  1  •  X . 


(2.14) 


Combining  equations  (2.12)  and  (2.14)  the  sum  of  the  forces 
in  the  x-direction  is 
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£F* 


3  3 

px  +  — (x  x  x  )  +  - ("^xy 

Sx  Sy 


dxdy . 


(2.15) 


The  rate  of  change  of  the  linear  momentum  in  the  x-direction 

assuming  constant  mass  is  mdu/dt. 

du  Su  dx  Su  dy  Su  Sx 

As  —  =  —  —  +  —  —  +  —  via  the  chain  rule,  and  —  =  u , 

dt  Sx  dt  Sy  dt  St  St 

Sy 

and  —  =  v  the  x-direction  change  in  linear  momentum  for 

St 

particle  is 


du  Su  Du  Su 

m —  =  pdxdy(v —  +  v —  +  — ).  (2.16) 

dt  dx  Sy  St 


Substituting  equations  (2.15)  and  (2.16)  into  the  x-component 
of  equation  (2.11)  yields 


Su 

pdxdy(u —  + 
Sx 


Su  Su 

v —  +  — ) 

Sy  St 


S 

px  +  — (txx) 

Sx 


S 

+  - (^xy) 

Sy 


dxdy . (2.17) 


Now,  in  order  to  have  the  entire  equation  as  a  function 
of  velocity  the  normal  and  shear  stresses  must  be  found  in 
terms  of  velocity. 

By  assuming  a  Newtonian  fluid  (Ref.  2:p.  B-5]  the  shear 
stress  is  linearly  related  to  the  rate  of  angular  deformation 
with  fluid  viscosity  being  the  proportionality  constant. 
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(2.18) 


Du  c)v 

"City  =  U( -  +  — — ) 

Dy  Dx 

The  normal  stresses  are  equal,  but  opposite  in  direction 
to  the  pressure  when  no  shear  stresses  are  involved.  With 
shear  stress  from  viscosity  it  is  assumed  that  the  normal 
stresses  deviate  from  -P  and  that  the  deviation  is 
proportional  to  both  a)  the  rate  of  linear  strain  in  the 
direction  concerned,  and  b)  the  rate  of  volume  deformation. 
Therefore,  the  normal  stress  in  the  x-direction  [Ref.  l:p.  B- 
10  ]  is 

Du  2  Du  Dv 

T**  =  -  P  +  2u(  — )  -  — Vl(—  +  — )  .  (2.19) 

Dx  3  Ox  Dy 

Applying  the  conservation  of  mass,  equation  (2.9),  equation 
(2.19)  simplifies  to 


Du 

=  -  P  +  2u(  — )  .  (2.20) 

Dx 


Substituting  equations  (2.18)  and  (2.20)  into  (2.17)  and 
dividing  by  dxdy  yields 


Du 

Du 

Du 

D  ciu 

D  ciu 

Dv 

p(u -  + 

V - 

+  __) 

=  pX  +  —  (-P  +  2 U(  — )  ) 

+  M — ( — 

+  — ) (2.21 

Dx 

Dy 

Dt 

Dx  Dx 

Dy  Dy 

Dx 

After  multiplying  and  rearranging  the  right  hand  side  becomes 
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Sp  S2U  c)2U  S2V 

px  -  —  +  2p -  +  P -  +  p - 

Sx  Sx2  c)y2  dySx 


which  is  also  equal  to 


Sp 

c)2U  S2u 

S  Su 

Sv 

px  - 

—  + 

p- —  +  P — 

+  p— ( —  + 

— ) 

Sx 

Sx2  Sy2 

Sx  Sx 

Sy 

Again  applying 

the 

conservation  of  mass 

,  equation 

(2.9)  , 

equation  (2.21) 

becomes 

Su 

Su 

Su 

Sp  S2u 

S2u 

p(u —  + 

V - 

+  _)  =  pX  - 

—  +  p(  — 

+  )  . 

(2.22) 

Sx 

Sy 

St 

ox  ox2 

Sy2 

With  v=  vi/p  =  kinematic  viscosity  and  neglecting  the  body 
force,  X,  the  two  dimensional  Navier-Stokes,  conservation  of 
momentum  equation  for  incompressible  and  constant  viscosity 
flow  in  the  x-direction  [Ref.  2:p.  436]  it 


Su 

Su 

Su 

i  Sp 

S2u 

S2u 

—  + 

u -  + 

V _  — 

—  _  _  + 

V  ( - 

+ 

St 

Sx 

Sy 

p  dx 

dx2 

Sy2 

Similarly,  in  the  y-direction  the  Navier-Stokes  equation  is 


Sv 

Sv 

Sv 

l  Sp 

S2v 

S  2  V 

—  + 

u -  + 

V -  = 

—  _  ....  + 

V  ( 

+  — ) 

St 

Sx 

Sy 

p  Sy 

Sx2 

Sy2 
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III.  STEADY  INVISCID  FLOW 


Although  real  fluids  are  viscous  the  major  effects  of 
viscosity  are  concentrated  in  a  region,  or  layer  adjacent  to 
a  body.  Therefore,  analyses  of  inviscid  flow  are  useful  and 
serve  as  a  good  approximation  to  flow  outside  the  boundary 
layer  and  wake  behind  the  body. 

The  justification  for  applying  the  results  of  perfect 
fluid  analyses  to  viscous  flows  was  postulated  by  Ludwig 
Prandtl  in  1904  [Ref.  3:p.  299].  He  stated  that  the  effects 
of  viscosity  on  the  flow  around  streamlined  bodies  at  high 
Reynolds  numbers  are  effectively  limited  to  a  "thin"  boundary 
layer.  The  characteristic  length  to  judge  thinness  is  the 
distance  from  the  forward  stagnation  point  to  the  point  of 
consideration . 

A.  VELOCITY  POTENTIAL 

For  flow  outside  the  boundary  layer  it  is  a  great 
advantage  to  simplify  equations  and  develop  a  single 
governing  equation.  With  the  assumptions  of  steady  flow,  no 
energy  transfer  to  or  from  the  fluid,  no  body  forces,  no 
shear  stress  (inviscid),  and  irrotational  flow  the  velocity 
potential,  4>/  is  utilized  [Ref.  3:p.  48].  4>,  a  scalar 
function  of  spatial  coordinates,  x  and  y,  is  defined  such 
that 

V  =  V4>  (3.1) 
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and 


0* 

c)4> 

a* 

c>y 

The  importance  of  the  velocity  potential  is  that  only  one 
equation  is  needed  to  describe  the  irrotational  flow. 
Velocity  components  are  obtained  using  equation  (3.2). 

B .  LAPLACE  EQUATION 

For  steady,  incompressible  flows  the  continuity  equation 
(2.9)  is 

c)u  dv 
c)x  c)y 

Rewriting  (2.9)  in  terms  of  the  velocity  potential  the 
equation  becomes 

dz4>  32<t> 

-  +  -  =  0  (3.3) 

c)x2  c)y2 

This  form  of  the  Laplace  equation  [Ref.  2:p.  81]  is  the 
governing  equation  for  steady,  irrotational  flow  of  an 
incompressible  fluid. 

The  importance  of  equation  (3.3)  is  that  it  is  linear 
allowing  for  the  principle  of  superposition.  For  example  if 
<K  /  4>z,  <^3 .  .  .  are  solutions  of  (3.3),  then  the  sum  c|>  =  4>x  + 

+  <t>3  +...is  also  a  solution  of  (3.3).  Superposition  of 
irrotational,  incompressible  flows  allow  for  the  construction 
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of  complex  flows  that  are  also  irrotational  and 
incompressible . 

C.  SIMULATION  (CONFORMAL  MAPPING) 

The  inviscid  flow  about  an  airfoil  can  be  obtained  most 
conveniently  by  means  of  a  transformation,  which  starts  with 
a  known  flow  about  a  simple  contour,  a  circle,  distorts  the 
contour  into  the  desired  shape,  and  simultaneously  adapts  the 
flow  to  that  shape.  The  transformation  is  accomplished  using 
a  sequence  of  three  conformal  mappings  [Ref.  l:p.  324]. 

The  first  mapping,  necessary  only  when  the  airfoil 
trailing  edge  has  non-zero  thickness,  is  accomplished  using  a 
logarithmic  mapping  function.  The  airfoil  is  perturbed 
slightly  to  make  the  upper  and  lower  surface,  trailing-edge 
points  coincide. 

The  second  mapping  analytically  removes  the  trailing-edge 
corner  using  the  Karman-Tref f tz  mapping. 

The  third  and  final  mapping  transforms  a  quasi-circular 
shape  into  a  perfect  circle  using  an  iterated  sequence  of 
Fast-Fourier  Transform  applications. 

During  the  transformation  of  streamlines  about  a  circle 
to  those  about  an  airfoil,  the  preferable  approach  insuring  a 
transformed  flow  free  from  vorticity  uses  the  complex 
variable  z  =  x  +iy  [Ref.  5:p.  285].  The  transformation  of  z 
to  another  plane  is  4=  f(z)  =  £  +  iT) .  The  potential 
function,  Q  =  4>  +  i'l',  is  irrotational  and  incompressible  in 
both  planes. 
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The  streamlines,  <l>  ,  and  equipotential  lines,  4>,  of  a  flow 
in  the  z-plane  will  transform  into  another  orthogonal  network 
of  lines  in  the  C -plane.  Different  magnification  ratios  and 
different  amounts  of  rotation  at  different  points  in  the 
field  will,  however,  change  the  appearance  of  the  flow 
pattern  from  about  a  circle  to  about  the  airfoil. 

The  general  transformation  function  whose  derivative 
dC/dz  satisfies  the  requirement  d£/dz  approaches  unity  for 
large  values  of  z  is 

C =  Z  +  ColnZ  +  Ca/Z  +  C2/z2  +...  .  (3.4) 

The  requirement  is  necessary  as  streamlines  are  not  distorted 
a  great  distance  from  the  body  where  the  body's  shape  has  no 
influence  on  the  flow. 

The  coefficients  may  be  real,  imaginary  or  complex.  A 
finite  number  of  the  coefficients  are  determined  from  the 
specified  normal  velocity  components  equally  spaced  around 
the  unit  circle,  and  from  the  Kutta  condition  which  ensures 
stagnation  at  the  trailing  edge. 

while  in  the  first  iteration  the  normal  velocities  are 
zero,  and  the  solution  for  flow  over  a  circle  is  used,  the 
subsequent  normal  velocity  boundary  conditions  are  determined 
from  the  previous  viscous-flow  calculations  using  the 
equation 
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d 

Vn  =  — (u„6~)  (3.5) 

ds 

where  u.  is  the  velocity  at  the  edge  of  the  boundary  layer 
and  6*  is  the  displacement  thickness.  Once  the  coefficients 
are  found,  the  real  and  imaginary  parts  of  equation  (3.4)  are 
equated  yielding 

K  =  £(x,y)  and  n  =  T)(x,y) . 

As  x2  +  y2  =  r2  the  two  equations  of  H,  and  r)  are  transformed 
to 


x  =  x(E,,r2)  and  y  =  y(t},r2). 


Then  x2  and  y2  are  added  to  yield 


x2  +  y2  =  r2  = 


=  x(  r2 )  2  +  y(ii,r2)  2.  (3.6) 


After  dividing  both  sides  by  r: 


1 

1 

r2 

x(4,r2) 

2  +  — 
r2 

Y(r\,r2) 

Then  each  circle  of  radius,  r,  in  the  z-plane  is  transformed 
to  the  proper  shape  in  the  £ -plane  to  describe  inviscid  flow. 
1.  Transformation  of  Velocities  [Ref.  5:p.  291] 

In  the  z-plane  as  Q(z)  =  4>(x,y)  +  i4»  ( x , y )  the 

velocities  are  defined  by 
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r 


dQ(z) 

-  =  V*  -  ivy  (3.8) 

dz 

where, 


V  =  V*  +  iVy. 


In  the  C -plane  the  velocities  are  also  found 


dQ 

—  =  V;  -  iVn  (3.9) 

dC  4 


where , 


vc  =  +  iVn  . 


The  velocities  in  the  two  planes  are  equated  by 


dQ  dz  V*  -  iVy 

V,  -  iV-  =  —  .  —  =  -  .  (3.10) 

1  dz  dC  dC/dz 


The  pressure  in  the  transformed  stream  is  related  to  the 
stream  velocity  through  Kelvin's  equation 


1  12 

-  p  Vr 1  2  +  p  =  constant  =  -p  Vr  +  pa.  (3.11) 

2  2^1 
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IV.  VISCOUS  FLOW 


A.  DERIVATION  OF  BOUNDARY  LAYER  EQUATIONS 

The  previous  analyses  provide  a  valid  solution  to  the 
flow  outside  the  boundary  layer.  Within  the  boundary  layer 


however. 

the 

effects  of 

viscosity  cannot 

be 

neglected. 

In 

laminar 

flow 

governing 

equations  can 

be 

obtained 

by 

simplifying  the  conservation  equations.  In  turbulent  flow, 
however,  the  number  of  variables  outnumbers  the  equations. 
Great  dependence  is  then  placed  on  dimensional  reasoning  and 
on  hypotheses  suggested  by  experimental  results. 

The  most  important  deduction  from  Prandtl's  thin  boundary 
layer  theory  is  that  static  pressure  can  be  considered 
constant  across  the  boundary  layer  [Ref.  3:p.  299]. 

Op 

—  =  o  (4.1) 

Oy 

As  the  boundary  layer  thickness,  6,  is  small,  d6/dx  is  also 
small.  Streamlines  are  then  only  slightly  curved  and  the 
radii  of  curvature,  R,  are  large.  With  a  large  R  the 
equilibrium  condition 

Op  u2 

Oy  P  R 
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illustrates  that  Dp/  Dy  will  be  very  small  and  can  be 
neglected.  Experimental  results  confirm  that  Dp/  Dy  may  be 
neglected  even  over  surfaces  of  small  radii  of  curvature. 

Also,  in  a  thin  boundary  layer  with  a  slowly  changing 
thickness  Dv/Dx  is  much  smaller  than  Du/Dy.  The  significance 
is  then  that  the  normal  shear  stress  may  be  neglected  when 
compared  with  the  viscous  shearing  stress.  Equation  (2.18) 
then  becomes 


Du 

Txy  =  P - .  (4.2) 

Dy 


With  this  simplification  the  approximate  equation  for 
flow  in  the  two-dimensional  boundary  layer  can  be  found 
directly.  Newton's  second  law  of  motion  applied  to  a  fluid 
element  of  mass  may  be  written 


Du 

Du 

Du 

Dp  Dx  y>« 

pdxdy ( — 

+  u - 

+  V— ) 

=  (-  —  +  - 

)  dxdy 

(4.3) 

Dt 

Dx 

Dy 

Dx  Dy 

as  illustrated  in 

Figure 

4.1. 

Substituting 

equation 

(4.2) 

and  dividing  both 

sides  by  dxdy 

yields 

Du 

Du 

Du 

Dp  D  Du 

p(—  + 

u —  + 

v— )  = 

(-  —  +  — (U— ) 

(4.4) 

Dt 

Dx 

Dy 

Dx  Dy  Dy 

In  terms  of  kinematic  viscosity  equation  (4.4)  becomes 


Du 

Du 

Du 

l  Dp  D2u 

—  + 

u -  + 

V -  = 

—  _  _ _ +  y  - 

Dt 

Dx 

Dy 

p  Dx  Dy2 
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Oyxdx 


Figure  4.1  Forces  Acting  on  an  Element  in  the  Boundary  Layer 
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This  equation  is  the  boundary  layer  equation  of  motion  and  is 
identical  to  the  equation  found  using  an  order-of-magnitude 
analysis  [Ref.  3:p.  443].  Equation  (4.5)  is  also  nearly 
identical  to  the  Navier-Stokes  equation  (2.23)  with  the 
exception  that  the  term  v32u/t)x2  is  deleted.  The  order-of- 
magnitude  analysis  suggests  that  this  term,  vc^u/Ox2,  may  be 
neglected  compared  to  vc)2u/  dy2 .  Combined  with  the 
continuity  equation  3u/c)x  +  dv/dy  =  0  (2.9),  equations  (4.5) 
and  (4.1)  are  known  as  Prandtl's  boundary  layer  equations 
[Ref.  2 : p .  P-9]. 

For  an  incompressible  flow,  there  are  three  variables,  u, 
v  and  p,  but  only  two  equations,  (4.5)  and  (2.9).  The 
equations  may  be  solved  though,  by  first  determining  p  as  a 
function  of  x  using  inviscid  methods,  setting  dp/c)y  =  0  in 
the  boundary  layer,  and  then  solving  (4.5)  and  (2.9)  for  the 
velocity  distributions. 

B .  TURBULENT  FLOW 

Turbulent  flow  as  differentiated  from  laminar  flow  is 
characterized  by  fluctuating  instantaneous  properties  which 
greatly  increase  the  complexity  of  the  problem.  A  very 
useful  simplification  to  the  turbulent  problem  is  then  the 
use  of  time-averaged  values,  denoted  by  a  bar  over  the  value. 
Instantaneous  values  are  indicated  by  the  prime  [Ref.  4:p. 
23  ] . 
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u  =  u  +  u’ 
v  =  V  +  V ' 
p  =  p  +  p' 

The  continuity  equation  containing  total  values  becomes 


3  _  0  _ 

— (u  +  U' )  +  - (V  +  V')  =  0. 

c )x  0y 


Simplifying  the  equation  becomes 


0  _  0  0  _  0 

— (u)  +  — (u' )  +  — (v)  +  — (v' )  =  0 

0x  0x  c)y  Sy 


with 


Z  _  0 

— (u)  and  — (u' ) 
c)x  0x 


0 

— —  ( u '  )  =  0 
0X 


0  _  0  _ 

—  (V)  =  - (V) 

0X  0X 


0 

and  — ( V  ) 
0X 


0 

=  - —  ( v  '  )  =  0. 
0X 


The  time-averaged  continuity  equation  for  turbulent  flow  is 
now 


0  _  0  _ 

— —  ( U )  +  —(v)  =  0  (4.6) 

0x  0y 
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Applying  total  values,  the  steady  version  of  the  Navier 
Stokes  equation  (2.23)  becomes 


_  3(u  +  u' )  _  c)(u  +  u' )  1  0(p  +  p' ) 

(u  +  u '  ) -  +  (v  +  v') -  =  - - 

Ox  Oy  p  Ox 


02(u  +  u' )  02 (u  +  u’ ) 


(4.7) 


After  simplifying,  the  time-averaged  Navier-Stokes  equation 
for  turbulent  flow  [Ref.  l:p.  C-10]  becomes 


1  Op  02u  02u 

p  Ox  Ox2  Oy2 


0  —  0  - 

-  — ~(u' 2)  -  — (U' V' ) .  (4.8) 

Ox  Oy 

The  new  terms,  0/Ox(u'2)  and  0/0y(u'v'),  which  correspond  to 
normal  and  shear  stress,  are  called  Reynolds 
stresses.  The  similar  y-component  terms  areO/Oy(v'2)  and 
0(v'u' )/0x. 

C.  TURBULENCE  MODELS 

The  time-averaged  Navier-Stokes  equation  is  nearly 
identical  to  the  original  equation  except  that  the 
instantaneous  values  are  replaced  by  the  mean  or  time- 
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averaged  values  and  two  additional  terms  involving 
fluctuating  velocities,  u'  and  v’ ,  appear.  An  interpretation 
of  these  two  terms  compares  them  to  the  previously  existing 
terms  c)2u/c)x2  and  Q2u/^y2.  The  right  hand  side  of  equation 
(4.8)  less  the  pressure  term,  and  after  multiplying  by 
density,  becomes 

c)2u  c)2u  £) -  - 

U- —  +  U -  -  P— u’  2  -  p  — u '  v ' 

c)x2  c)y2  c)x  c)y 


or 

c)  c)u  -  c)  c)u  - 

—  (U—  -  pu'2)  +  — (p —  -  p  u  ’  v '  )  . 

Ox  c)x  <3y  c)y 

As  each  term  has  the  dimensions  of  stress,  and  p(c)u/c)y)  is 
part  of  the  laminar  shear  stress  xv«,  it  appears  that  the 
term  -pu'v'  represents  a  turbulent  addition  to  shear  stress 
[Ref.  2 : p .  T-2 ] .  Now,  this  shear  stress  is  really  a  vertical 
mixing  of  horizontally,  travelling  fluid  particles.  A  model 
of  this  mixing  then  calculates  the  rate  of  momentum  transfer 
involved. 

1.  Prandtl's  Mixing-Length  Model 

To  predict  the  turbulent  stresses  Prandtl  assumed 
that  turbulent  fluctuations  are  primarily  the  result  of  the 
mean  velocity  differences  between  two  layers  in  the  flow. 
Therefore,  u’  is  proportional  to  Ou/c)y  with  the 
proportionality  factor  having  the  unit  of  length. 
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(4.9) 


-  3U 

(u12)*  =  a — 

By 


Also,  assuming  that  v'  is  of  the  same  order  of  magnitude  as 
u'  at  a  particular  point, 


-  Bu 

(v'2)*  =  b — .  (4.10) 

By 


Substituting  u'  and  v’  into  the  turbulent  shear  stress,  x-r, 
is 


-  Bu  Bu 

-  pu'v’  =  -  pab( — )( — ).  (4.11) 

By  By 

As  a  and  b  are  both  unknown  constants  of  length,  they  both 
may  be  replaced  by  the  "mixing  length",  1,  the  hypothetical 
distance  between  the  two  layers  involved.  The  turbulent 
shear  stress  is  now 


T«r 


Pi2 


Bu  Bu 

By  By 


(4.12) 


which  insures  the  correct  sign. 

2.  Eddy  Viscosity  -  Cebeci-Smith 

The  turbulent  addition  to  shear  stress  may  also  be 
modeled  in  terms  of  "eddy  viscosity".  As  laminar  s’ ear 


25 


stress  =  Ty><  =  u(du/c)y)  =  pv(c)u/  dy,  turbulent  shear  stress 
may  be  equated 


-  c)u 

-  p u 1 v '  =  pvt — ,  (4.13) 

Qy 

where  vt,  the  eddy  viscosity  is  empirically  determined. 

The  method  used  here  is  that  of  Cebeci  and  Smith  as 
modified  by  Cebeci,  Clark,  Chang,  Halsey  and  Lee,  [Ref.  l:p. 
327] . 

Eddy  viscosity  by  Cebeci  is  represented  by 


{0 . 4y[ 1  -  exp  (-  -) ]}2 

A 


Ttr  (0  <  y  <  yc )  (4.14) 


a 


(ue  -  u)dy 


Ytr-Y  (yc  <  y  <  6)  (4.15) 


where 


A 


26v 


and 


1 

Y  =  -  . 

1  +  5 . 5 (y/6) * 


The  distance  from  the  body  yc  which  is  less  than  the  boundary 
layer  thickness,  6,  is  the  distance  where  the  two  equations 
(4.14  and  (4.15  give  the  same  resultant  vt. 
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The  intermittency  factor,  Ytr.,  which  indicates  the 
local  fraction  of  turbulent  flow  to  total  flow,  is  given  by 


The  location  of  the  start  of  transition  is  xtr,  and  the 
empirical  factor  G  is 

1  u„3 

G  =  -  -  R-1 ' 3*  (4.17) 

1200  v2 

where  Rxtr  is  the  transition  Reynolds  number 

Rxtr  ~  ( UeX/v  ) tr  •  (4.18) 

In  equation  (4.17)  the  empirical  constant  GYtr-  = 
1200,  was  used  by  Chen  and  Thyson  [Ref.  4:p.  327],  Values 
lower  than  1200  may  give  better  results  at  low  Reynolds 
numbers  as  will  be  discussed  in  Section  VII. 

The  expression  for  a  is 

0.0168 
a  =  - 

p2  .  S 

The  non-dimensional  factor  F  represents  the  ratio  of 
the  product  of  the  turbulent  energy  by  normal  stresses  to 
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that  by  shear  stress  evaluated  at  the  location  where  shear 
stress  is  maximum. 


F 


V  2)c)u/c)x 


-  u'v'  c)u/c)y 


( -  U'V') max 


(4.20) 


According  to  the  data  of  Nakayama  [Ref.  l:p.  327],  £  can  be 
represented  by 


£ 


-  u'v' 


(  -  U'V') max 


6 

. .  11  R*r  ^  1 

1  +  2Rt(  2  *  R*r ) 


1  +  Rt 

-  Rx  >  1 

Rt 


where  RT  =  xw/ ( -  u'v')Max  and  xw  is  the  wall  (body)  shear 
stress . 


D .  TRANSITION 

The  location  of  the  onset  of  laminar-to-turbulent 
transition  when  not  found  experimentally  is  determined 
empirically.  The  method  used  by  Cebeci  [Ref.  4:p.  333]  is 
the  criterion  proposed  by  Michel. 

At  the  point  of  transition  the  Reynolds  number  based  on 
momentum  thickness,  ©,  is  related  to  the  Reynolds  number 
based  on  the  coordinate  position,  x. 
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22400 

Retr  =  1.174  (1  +  - ) R.xtr  °*6  (4.21) 

R»ictr 


where 


Hex  —  Ue(X/v )  and  Retr  —  Ua(0/v) • 


In  the  Cebeci  Code  the  transition  may  be  determined  in 
the  following  ways: 

1)  The  points  of  transition  are  calculated  using  Michel's 
criterion . 

2)  If  laminar  separation  occurs  forward  of  the  criterion 
points,  Michel's  criterion  is  disregarded  and 
transition  is  redefined  at  the  separation  point. 

3)  The  transition  locations  may  be  specified  by  the  user  . 
provided  stall  is  not  computed. 
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V.  VISCOUS  METHODS 


Momentum  transfer  in  fluids  is  accomplished  by 
hydrostatic  pressure  and  viscous  stresses.  When  viscous 
stresses  are  negligible,  fluid  behavior  can  be  predicted  by 
inviscid  flow  methods  as  stated  in  Section  III. 

Viscous  stresses  caused  by  a  variation  in  velocity  in  a 
direction  normal  to  the  flow  are  called  shear.  The  most 
common  shear  is  that  found  in  the  boundary  layer  between  a 
displayed  stream  and  the  solid  surface.  With  the  "no  slip" 
condition  fluid  velocity  is  zero  on  the  surface,  but  the 
velocity  gradient  is  not  so  constrained.  From  the  body  along 
its  normal  direction  the  fluid  velocity  asymptotically 
approaches  that  of  the  free  stream. 

As  mentioned  in  Section  III,  Prandtl  hypothesized  the 
division  of  the  flowfield  into  the  two  regions,  the  boundary 
layer  where  viscous  effects  cannot  be  neglected  and  the 
region  outside  the  boundary  layer  which  may  be  considered 
inviscid. 

This  hypothesis  allows  for  the  use  of  the  parabolic 
boundary  layer  equations  of  section  III  instead  of  the 
elliptic  Navier-Stokes  equations.  Depending  on  the  boundary 
conditions,  solutions  fall  into  three  methods  [Ref.  6:p.  13]: 

1)  The  direct  boundary  layer  method.  This  method  uses  the 
"no  slip"  condition,  where  normal  and  tangential 

velocities  are  zero  at  the  surface,  and  a 

pre-determined  velocity  is  specified  at  the  boundary 
layer  edge. 
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2)  The  inverse  boundary  layer  method.  Boundary  conditions 
are  replaced  by  wall  shear  or  displacement  thicknesses. 

3)  The  interactive  boundary  layer  method.  The  edge 
boundary  condition  drives  a  combination  of  displacement 
thickness  and  external  velocity. 

Methods  one  and  three  will  be  discussed  as  they  apply  to 
the  Cebeci  Code. 

A.  DIRECT  BOUNDARY  LAYER  METHOD  [Ref.  6:p.  13] 

This  method  for  solving  the  boundary  layer  equations  is 
used  only  near  the  leading  edge  where  the  viscous  effects  are 
small.  Initial  conditions  are  generated  at  the  stagnation 
point  and  the  equations  are  integrated  around  the  leading 
edge.  The  numerical  solution  utilizes  a  finite  difference 
method  where  the  continuity  and  momentum  equations  are 
redefined  as  a  system  of  linear  algebraic  equations. 

The  method  begins  by  describing  steady,  incompressible, 
2-D  flows  in  a  curvilinear  coordinate  system  where  x  is 
directed  along  the  airfoil  surface  and  y  is  perpendicular  to 
x.  The  boundary  layer  equations  with  the  turbulent  Reynolds 
stress  are 


du  dv 

—  +  —  =  o  (5.1) 

dx  dy 


1  dp  l  d  du  -  _du  _du 

—  + - [y —  -  pu’v'  ]  =  u —  +  v —  (5.2) 

P  dx  p  dy  dy  dx  dy 
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0. 


(5.3) 


dp 
dy 

where  the  order  of  magnitude  of  any  turbulent  stress  is 
assumed  to  be  that  of  its  laminar  stress.  The  boundary 
conditions  are: 

aty=0;u=0,v=0  (5.4) 

at  y  =  « ;  u  =  U* ( x )  (5.5) 

The  eddy  viscous  stress  previously  defined  is  reprinted  as 

-  du 

-pu'v'  =  pvt -  (4.13) 

dy 

Also,  the  pressure  gradient  term  may  be  written 

1  dp  due 

--  —  =  u»_  (5.6) 

p  dx  dx 

Therefore,  the  momentum  equation  (5.2)  may  be  rewritten  as 
_du  +_  du  du«  d  du 

u—  V— —  =  u. -  +  —  (b— )  (5.7) 

dx  dy  dx  dx  dy 

where  b  =  1  +  v*/v  and  the  boundary  conditions  are 

at  y  =  0;  u(x,0)  =  0,  v(x,0)  =  0  (5.8) 

at  y  =  y«;  u(x,y.)  =  u«(x) . 
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1.  Falkner-Skan  Transformation  [Ref.  6:p.  14] 


To  solve  the  new  boundary  layer  equations  the 
Falkner-Skan  transformations  are  used,  which  reduce  the 
number  of  variables,  and  scale  the  normal  coordinate  y  and 
the  stream  function  4'  with  reference  to  the  external 
velocity. 


T)  = 


(5.9) 


4»(x, y)  =  ( uB  x)*  •  f  ( x ,r)  )  (5.10) 

The  continuity  equation  is  automatically  satisfied 
using  the  stream  function  (u  =  d|/0y  and  v  =  -O’l'/Ox)  . 
Therefore,  only  the  momentum  equation  needs  to  be  solved, 
which  after  transformation  becomes 

m  +  1  Of '  Of 

(bf")’  +  -  ff"  +  m [ 1  -(f)2]  =  x ( f ' —  -  f" — ) ( 5 . 1 1 ) 

2  Ox  Ox 

where  m  =  ( x/ue ) ( duo/0x) ,  a  dimensionless  pressure-gradient, 
and  f'  =  Of /Orj. 

This  equation  (5.11)  is  a  third  order  partial 
differential  equation,  and  the  solutions  are  ”non-similar"  as 
they  are  functions  of  both  x  and  i) .  If  the  solutions  were 
only  a  function  of  i)  ,  then  the  right  hand  side  of  the 
equation  would  equal  zero,  and  the  flow  would  be  "similar" 


33 


[Ref.  1 : p .  v-10].  To  solve  equation  (5.11),  a  numerical 

solution  is  needed  such  as  the  "box"  method. 

2 .  The  Box  Method 

The  box  method,  developed  by  Keller  in  1970  [Ref. 
l:p.  331],  is  a  widely  used  methods  for  solving  non-linear 
differential  equations.  The  steps  of  this  method  include  the 
conversion  of  the  Falkner-Skan,  transformed,  momentum 
equation  into  a  system  of  first-order  partial  differential 
equations.  This  non-linear  system,  after  conversion  from  a 
continuous  function  to  discrete,  is  linearized  by  Newton's 
method.  The  block  elimination  method  is  then  used  to  solve 
the  linearized  difference  equations  of  the  boundary  layer 
problem . 

The  third  order  momentum  equation  (5.11)  is  converted 
into  a  first  order  system  with  the  addition  of  the  dependent 
variables  U  and  V  [Ref.  6:p.  14]. 

U  =  f'  (5.12) 

V  =  U'  =  f"  (5.13) 

m  +  1  Ou  Of 

(bV) '  +  - fv  +  m( 1  -  U2)  =  x(U —  -  V — )  (5.14) 

2  Ox  Ox 

The  boundary  conditions  are 

at  T)=  0;  U(x,0)  =  0,  f ( x, 0 )  =  0 
at  i)  =T}«;  U(x,t}*)  =  1 
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The  solution  domain,  0<x<xT  and  0<n*T)a,  of  the 
continuous  functions  f,  U  and  V  is  then  covered  by  a 
rectangular  grid  to  facilitate  the  problem  solving  with  a  set 
of  discrete  values.  This  grid  is  shown  in  Figure  5.1. 

The  subsequent  notation  [  is  used  to  represent  the 
quantities  of  f,  U  or  V  at  the  point  (x^,  1)3  )■ 

All  quantities  can  then  be  approximated  by  network 
coordinate  values.  Using  the  grid  system,  the  solution  of 
the  parabolic  layer  equation  at  a  certain  streamline  position 


depends  solely  on  the 

solution  of 

upstream 

positions , 

while 

no  downstream  influence 

needs 

to  be 

considered. 

As 

calculations  proceed 

from 

the 

stagnation 

point  in 

the 

downstream  direction,  the  overall  solution  can  be  obtained 
incrementally.  Hence,  one  step  of  the  solution  procedure 
sets  up  the  governing  equations  for  a  column  of  grid  boxes  in 
the  sub-domain 


Xi.j<x<Xi  and  0<ti<tij 

and  solves  for  the  values  of  the  downstream  grid  position. 
The  x-grid  position  currently  solved  for  is  then  assigned  the 
superscript  "i"  while  "i-l"  represents  the  previous  position 
of  known  properties.  Using  coordinates  of  box  midpoints  and 
centered-dif f erence  derivatives,  the  equations  are  actually 
satisfied  midway  between  the  grids. 

Equations  (5.12)  through  (5.14)  in  terms  of  finite 
differences  [Ref.  6:p.  15]  are  now  written 
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LEGEND 

O  known 

□  unknown 

A  center  for  momentum 
equation 

V  center  for  equations 
containinq 
ffderivatives  only 


Figure  5.1  Rectangular  Grid  for  Finite  Difference 

Approximation 


36 


(5.15) 


+  U1 
D 


U3  -  U3-l 
- 


?  <v5  +  vj-i> 


V)*"1  -  (bV)j'1 


l- 

j-i 


l- 

j-J 


(5.16) 


(5.17) 


where  the  ordinary  differential  equations,  (5.15)  and  (5.16), 
are  centered  about  (xi7  nj_i)  and  the  partial  differential 
equation,  (5.17),  is  centered  about  (Xi_4,  r^-*). 

The  boundary  conditions  are 

U*  =  0,  f*  =  0  and  uj  =  1 

Equations  (5.15)  and  (5.16)  are  the  centered  difference 
derivatives . 

3.  Newton's  Method  [Ref.  6:p.  15] 

This  set  of  finite  difference  equations  is  nonlinear 
with  combinations  of  unknowns.  Newton's  method  is  therefore 
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needed  to  solve  the  system.  The  variables  are  linearized 
using  an  iterative  procedure  with  preceding  values. 

fj'1  =  fj”1+  6fj  1  and  fj,K  =  fj'K_1  +  6fj‘K  for  Kl2 

i  r  r  1  •  K  .  K-l 

where,  6f^  << 

U* ' 1  =  U^"1+  6U j ‘ 1  and  U*'K  =  Uj'K_1  +  6^‘K  for  K>2 
where,  6uJ,K<<  U^-*"1 

Vj  1  =  vj”1+  ^Vj'1  and  V*'K  =  V*,K_1  +  6V j ' K  for  K>2 
where,  6Vj'K<<  Vj'K  1 

K  is  the  iteration  counter.  After  substituting  these  values 
into  equations  (5.15)  through  (5.17)  and  neglecting  the 

i , K  i , K  i.K 

quadratic  terms  of  6f  ,  6U  and  6V  ,  the  system  of 

j  j  j 

unknowns  is  then  linear. 


(5.18) 


(5.19) 
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(SO 


t  a  \  X  •  X  *  x\ 

+  (s2)  .  eVj 


+  + 


(so  * ,Ksf *1*  +  (SO  j-K6uj-K  +  (S6)j-K6uj;^ 


(r2)j‘K 


C 


where 


hi.K-l  x 

(Sx)^K=  jh  —  +  ^  (f^.'K“1 

J  A*-i 


fj-r  3-1 


mi+1  O-K-l 


'  j-i 


(Sa)  J'K= 


bi-K-l  x 

bj  xi-|  (fi.K-l 

TT~  (fj-| 


3 


-i-1.  ,  m  +  1  -i.K-l 

fj-r  t~ 


3-1 


(S3)i-K=  ^  (vili"1  +  v;_t)  + 


j  -  nq-  'wj- 


i-1,  .  m1  +  1  ..i.K-1 

3  Vj 


3-i 


/  q  <  i  •  3  -  5  mO  ’  K  1  A  \J  1  ) 

{S*]i  -  TR—  (VjO  V3-r 


nr  +  1  ,.i  .  K-l 
+  - 5—  V3-l 


(s5)  j,K= 


x . 


-  (  +  mi)U^K‘1 

Ki  3 


(so  yK- 


x. 


,  A  i .  i . K-l 

<-0  *  “  )uj-i 


,  ,i-K 

( ra ) j  = 


/  U.»»  \  i  •  K-  1  /  V\17  \  i  •  K"  1 

(bV)  .  -  m  +  1  i  K-l 

<( - 3 - j_ - Li_  ♦ 


,xi-i  *  i /  r.i •  K-l , z  ^  *1-1, „1. K-l  fi.K-l 


X, 


-  vj:j  fj:r 


-  {(■ 


(bv) 


l-f 


D  ^  i  —  }  * 

-  (bV)3-I  .  mi-1  +  1 

TT -  +  — -2 - 


fV) 


i-1 

3-\ 


+  m 


i-1 


(U 


i-1 

3-i 


(V 


i-1 

j-4 


f 


i-1 

j-i 


+  2m1"*} 


The  boundary  conditions  are 

6f o  =  0,  6U0  =  0,  6U^  =  0 
4.  Block  Elimination  Method  [Ref.  6:P.  17] 

The  system  of  equations  are  iteratively  solved  until 
i.K  i.K  i.K 

6f  ,  6U  and  6v  become  small  enough  to  be  neglected. 

3  3  3 

The  solution  method  is  that  by  Keller  and  is  called  the  Block 
Elimination  Method.  In  this  method  block-tridiagonal 

matrices  are  composed  of  blocks.  Only  those  blocks  on  the 
main  diagonal  and  on  the  two  ad3acent  diagonals  are  non-zero. 
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[aJ-k]  [cj'K] 


1 

i.K 


[b*'*]  •  k]  [c| * k] 


[B j " K ]  [A*'*]  [Cj-K] 


{6^*K> 

{6j'K> 


{6 j ' K} 


ir['K) 

{r^K} 


{r j  *  K> ( 5 . 21 ) 


rBi-K 

[BJ-1 

]  [aJ;*] 

icj:*) 

/  c  i  •  K  •> 
{6J-1} 

— 

[Bj*K] 

[*j'k£ 

{6j ‘ K> 

where  the 

blocks  are 

3x3  matrices 

T 

-hj/2 

0 

— 

[a^-k]  = 

( s3 

,,K 

(Ss)  j’K 

(Si)  j 

.K 

2< j  < J-l 

0 

-1 

! 

“1 

-hj/2 

0 

_ 

[b^-k]  = 

(S4 

»}-K 

(S.)  j’K 

(Sz)  j 

.K 

2<j<J 

0 

0 

0 

15 

0 

0 

— 

[cj'K]  = 

0 

0 

0 

0 

1 

"hj+l 

/2 

T 

0 

0 

— 

[  A  j  '  K  ]  = 

0 

1 

0 

0 

-1 

-h2/2 
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[A 


1  -hj/2  0 

( S  3 ) j " K  (Ss)j,K  ( S  a  ) j  *  K 

0  10 


The  right-hand  sides  of  the  equations  are 


.  .  i . K  _  -i.K-1  fi . K-l  .  TTi  •  K-l 

(Fl)j  '  fj-l  '  fj  +  hjUj-l 


(r2)j‘K  =  {As  defined  in  equation  (5.20)} 


(r3) 


i.K  „i . K-l  „i.K-l 


u 

3 


i.K-1 


uj-l  +  hj+lvj+i 


2<h<  J 


2<3<J 


2<h<  J 


(r)^,K  =  0,  ( r  2 ) ^ ’ K  =  0-  =  0 


The  block  elimination  method  solves  the  linear 
equations  with  two  steps.  The  forward  step  eliminates  the 
lower  diagonal  of  the  tridiagonal  matrix.  The  reverse  step 
solves  the  remaining  system  from  bottom  to  top. 

B.  INTERACTIVE  BOUNDARY  LAYER  METHOD  [Ref.  6:p.  18] 

As  the  direct  boundary  layer  method  is,  as  previously 
stated,  restricted  to  regions  of  small  viscous  effects,  and 
integration  of  the  boundary  layer  equations  fails  at  points 
of  zero  skin  friction,  a  method  is  needed  to  integrate  the 
boundary  layer  through  the  point  of  emerging  reversed  flow. 
This  method  must  also  account  for  strong  interaction  effects 
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due  to  separation  and  rapid  acceleration  of  the  flow 

downstream  of  the  trailing  edge. 

The  interactive  method  fulfills  these  requirements  by 
treating  the  external  velocity  and  displacement  thickness  as 
unknown  quantities.  Reflecting  the  elliptic  nature  of  the 
outer  flows,  an  additional  unknown  is  introduced,  but  the 
solution  can  be  obtained  using  either  the  eigenvalue  or 
Mechul  function  methods. 

The  Mechul  method  is  preferred  as  the  eigenvalue  method 
involves  nonlinear  problems.  In  this  method  the  edge 

boundary  condition  of  the  direct  method  is  supplemented  with 
the  interactive  boundary  condition.  The  unknown  external 
velocity  is  related  to  its  displaced  and  perturbed 
conditions.  The  unknown  functions  u(x,y),  v(x,y)  and  ue(x,y) 
are  represented  in  this  system  of  boundary  layer  equations 

c)u  3v 

—  +  —  =0  (5.22) 

dx  c)y 
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Que 
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0y 

where  pressure  in  the  y-momentum  equation  is  expressed  in 
terms  of  the  external  velocity. 

The  Mechul  function  approach  assumes  that  the  external 
velocity,  ue,  is  a  function  of  two  arguments,  x  and  y, 
allowing  for  an  easy  setup  of  finite  difference  equations, 
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and  avoiding  nonlinear  eigenvalue  techniques.  The  velocity 
components  u  and  v  are  required  to  satisfy  the  no-slip 
condition  at  the  surface,  and  u  must  merge  smoothly  into  the 
edge  velocity. 

at  y  =  0:  u(x,0)  =  0,  v{x,0)  =  0 

at  y  =  y„:  u(x,Y»)  =  u„(x,Y«),  u0(x,ye)  = 


=  Uel 


d  dE. 

— ( ue6* ) - 

d£  x-4 


Ue x ( x )  is  the  inviscid  edge  velocity  and  the  last  term, 
the  Hilbert  integral,  approximates  the  viscosity  induced, 
perturbation  velocity. 

Interactive  methods  are  useful  in  both  attached  and 
separated  regions,  while  direct  methods  fail  at  the  onset  of 
reversed  flow,  and  inverse  methods  converge  poorly.  Only  at 
the  stagnation  point  singularity  r^e  interactive  methods 
prohibited.  The  transformation  of  the  partial  differential 
equations  into  a  linear  system  of  algebraic  equations  is  very 
similar  to  that  of  the  direct  method.  The  normal  coordinate 
y,  streamfunction  <J< ,  and  the  external  velocity  u„  are  scaled 
with  reference  to  a  constant  velocity  u0,  and  the  local 
streamwise  coordinate  x. 


n  = 


(5.25) 
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4»(x#y)  =  (u0vx)  *  f  (x,r)  ) 


(5.26) 


ue (x, y) 

W(x,ij)  =  -  (5.27) 

Uo 

U0  is  the  vector  mean  velocity.  Ue  cannot  be  the  reference 
velocity  as  for  Falkner-Skan  variables  because  in  this  case 
the  external  velocity  is  unknown.  The  first  order,  semi- 
transformed  coordinate  system  with  additional  variables  U  and 
V  is 


U  =  f' 


(5.28) 


V  =  U'  =  f"  (5.29) 


1  c)w  c)u  c)f 

(bv)'  +  -  fv  +  XW -  =  X ( U -  -  V - )  (5.30) 

2  dx  <Dx  c)x 


W'  =  0 


(5.31) 


and  the  boundary  conditions  are 


at  f)  =  0 : 
at  T)  =  r)e 

W(  X,  T)«  )  — 


U(  x ,  0  ) 

:  U(x,tu) 

ueI(x) 

-  + 

Uo 


=  0,  f(x,0)  =  0 
=  W(X,1J.) 


—  {(^)*(  W(^ ,  n-)  n- 

d£  Uo 


-  fU/T).)]}- 


x-i. 
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The  conversion  of  the  flowfield  to  discrete  values  is 
very  similar  to  that  of  the  direct  method  with  an  orthogonal 
grid,  central  differences,  and  two-point  averages.  In  the 
interactive  method  though,  the  solution  proceeds  in  the 
downstream  direction  only.  As  only  downstream  disturbances 
are  accounted  for,  backflow  causes  numerical  instabilities. 
A  stable  integration  can  be  accomplished  though,  with  the 
assumption  that  backflow  velocities  are  comparatively  small. 
The  FLARE  approximation  (Flugge-Lotz  and  Reyhner) ,  [Ref.  6:p. 
19]  sets  the  streamwise  convection  term  ucWc)x  equal  to  zero 
in  regions  of  backflow. 

1  if  U1  i  >0 
D“2 

0  if  u1  !  <0 

3~  2 

The  finite  difference  equations  of  the  interactive 
boundary  layer  with  the  "on-off  switch"  are  now 


FLR 
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3  3 


(5.32) 
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=  X 


i-I 


fLR3- 


i-i  u1  -  u1”1 

I(U  K  '  j-i 


-  (V 


i-i  f' 


-  f 

■*r 


i-i 


,1 


0. 


(5.34)  and  (5.35) 


The  boundary  conditions  are  also  expressed  in  terms  of 
grid  or  nodal  values.  A  panel  method  type  approximation 
leads  to 


uj  =  0,  fj  =  0 


9i  +  cii<wj"j 


where  g±  and  ctl  represent  a  parameter  and  the  diagonal 
element  of  the  interaction  matrix,  due  to  a  discrete 
approximation  to  the  Hilbert  integral.  To  keep  the  number  of 
generated  terms  to  a  minimum,  ordinary  differential  equations 
like  the  y-momentuin  equation  are  centered  about  the 
downstream  face,  and  partial  differential  equations  like  the 
x-momentum  equation  are  centered  about  the  middle  of  the  box. 

The  unknowns  occur  in  vectors  of  four  components 


{fj'  Uj  '  Vj'  Wj}‘ 

The  J  quadruplets  of  unknowns  match  with  4J  equations, 
including  2(J-1)  auxiliary  relations  and  (j-1)  x-momentum  and 


(j-1)  y-momentum  equations.  Each  equation  corresponds  to  one 
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of  the  ( j  —  1 )  grid  rectangles  and  four  boundary  conditions. 
The  system  is  linearized  around  the  values  of  the  preceding 
iteration  (counter  K-l)  and  a  system  with  Newton  iterates 


c*i.K  c„i.K  Ctri.K  c„i.K 

6fj  /  ,  6Vj  ,  6Wj  ,  emerges. 


6f*'K  - 


r2(6uJ-K  +  6a*;*) 


f i . K-l  f i . K-l 

j-1  "  j 


+  h3  “i-f1 


(5.36) 


6U  ^  '  K  -  6U*;* 


-^i(6V^‘K  +  6vV*) 


1  -  u*-**1 


+  h3 


(5.37) 


(Si)J,K6Vj,K  +  (S2)j’K6V*;*  +  (S3)J,K6fJ'K  + 


(S„)*-K6f*;*  +  ( Ss )  j  ‘  K6U  ^  ‘  K  +  (S6)JlK6U*;*  + 


(  S7  )  j  ’  K6W  j  '  K  +  (sB)j'K6wJ;*  =  6wJ;*  =  (rz)j*K 


(5.38) 


6Wi-K 
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(5.39) 
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The  boundry  conditions  are 

6uJ-K  =  0,  6f^*K 


=  0,  6U j  * K  -  6Wj‘K  =  0 


cfi.K  x  ,1  „  ...i.K  yi  .i.K-1  .1  „  .  „i.K-l 

j  cii  J  J  cii  J  cii  J  J 


The  overall  procedure  is  a  repetitive,  linear  approach  to 
solve  the  nonlinear  system.  The  numerical  solution  is  again 
obtained  using  the  block  elimination  method  by  Keller,  except 
that  unlike  the  direct  method  the  vectors  of  the  unknown 
Newton  iterates  are  four  dimensional. 


<6£-K( 
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In  matrix-vector  form 
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The  4x4  matrix  blocks  are 
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(r2)^'K  =  {As  defined  in  equation  (5.38}} 
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VI.  INTERACTION  METHODS  [Ref.  7:p.  79 J 


Interactive  methods  couple  viscous  and  inviscid  flows  and 
are  intended  to  compute  through  regions  of  flow  separation. 
Given  their  levels  of  success,  these  methods  have  become 
inexpensive  alternatives  to  the  Navier-Stokes  solvers. 

The  simple,  classical  method  of  computing  the  viscous 
flow  over  an  airfoil  is: 

1)  The  velocity  distribution  is  computed  for  inviscid 
flow. 

2)  The  inviscid  velocity  is  input  to  the  viscous  flow. 

3)  The  viscous  flow  is  computed  by  integrating  the 
boundary  layer  equations. 

Now,  this  method  is  good  at  predicting  lift  and  drag,  but 
only  if  the  flow  remains  attached,  as  information  is 
transferred  only  once  from  inviscid  to  viscous  regions.  For 
more  complex  flow  multiple  information  transfers  are 
required . 

Close  coupling  is  needed  to  compute  flows  with  separation 
or  separation  bubbles.  A  better  method  than  the  previously 
outlined  classical  method  for  exchanging  information  between 
viscous  and  inviscid  regions  is  interaction.  The  different 
elements  of  interaction  include  direct  and  inverse,  inviscid 
and  viscous  flow  solvers.  Table  6.1  illustrates  the 
different  elements. 

The  disadvantage  of  the  direct  boundary  layer  method  is 
that  the  equations  become  singular  at  the  point  of 
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separation.  The  point  of  separation  may  be  integrated 
through,  however,  if  the  external  velocity  is  computed  with  a 
predetermined  displacement  thickness.  This  method  is  known 
as  the  inverse  boundary  layer  method. 

Another  problem  associated  with  separation  is  the 
instability  of  numerical  methods  which  prohibits  downstream 
marching  in  regions  of  reversed  flow.  In  the  situation  where 
flow  is  reversed  the  FLARE  approximation  is  used,  where  the 
momentum  transport  term  u3u/c)x  is  neglected.  This 
approximation  is  not  necessarily  accurate,  but  it  does  allow 
for  continued  calculations. 

Four  interaction  models  have  been  developed  to  calculate 
combined  inviscid  and  viscous  flows.  All  procedures  solve 
the  Laplace  equation  for  inviscid  flow  and  the  boundary  layer 
equations  for  viscous  flow.  The  four  models  are  the  direct, 
inverse,  semi-inverse  and  viscous-inviscid  interaction 
methods.  Each  model  is  subject  to  different  boundary 
conditions . 

The  first  three  models  are  considered  weak  interaction 
methods  in  that  they  provide  only  a  loose  coupling  between 
viscous  and  inviscid  regions.  The  two  regions  are  treated 
alternately.  As  indicated  in  Table  6.1,  the  viscous  flow 
solver  calculates  the  flow  in  the  viscous  region  and  produces 
the  boundary  condition  of  the  inviscid  region.  The  inviscid 
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TABLE  6 . 1  INTERACTION  ELEMENTS 


BOUNDARY  ( 

CONDITION 

Flow 

Direct 

Inverse 

Inviscld 

*  Zero  normal 

*  Prescription  of 

velocity  at 
the  surface 

velocity  distribution 

Viscous 

*  No  slip 
condition 

*  No  slip  condition 

*  Prescription  of 

*  Prescription  of  dis- 

external  velocity 

placement  thickness 

flow  solver  calculates  the  flow  in  the  inviscid  region  and 
produces  the  boundary  condition  of  the  viscous  region.  The 
weak  interaction  methods  process  either  displacement 
thickness  or  external  velocity  as  an  input  and  the  other 
quantity  as  an  output. 

In  contrast,  the  fourth  method,  the  simultaneous 
interaction  method,  is  considered  a  strong  interaction 
method.  A  strong  method  calculates  displacement  thickness 
and  external  velocity  simultaneously.  The  foundation  of  the 
four  interaction  methods  are  discussed  below. 

A.  DIRECT  INTERACTION  METHOD 

The  direct  interaction  model  is  composed  of  direct 
inviscid  and  viscous  flow  solvers  as  indicated  in  Figure 
6.1a.  The  external  velocity  distribution  is  calculated  first 
by  inviscid  computations.  The  displacement  thickness,  6*,  is 
then  calculated  using  the  external  velocity  as  a  boundary 
condition.  An  updated  shape  of  the  displacement  body  is  then 
computed,  and  all  steps  are  recomputed  in  order  until  the 
results  converge.  As  previously  stated  this  method  breaks 
down  at  the  point  of  separation,  and  is  therefore  not  useable 
in  regions  of  separated  flow.  However,  it  is  very  useful 
where  viscous  effects  are  small.  The  direct  method  is  used 
in  the  Cebeci  Code  around  the  nose  and  stagnation  point  of 
the  airfoil. 
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B.  INVERSE  INTERACTION  METHOD 


T?'is  method  was  developed  to  circumvent  the  singularity 
problems  near  separation.  According  to  Figure  6.1b  it  uses 
inverse  inviscid  and  viscous  flow  solvers.  Because  of  the 
inverse  method’s  very  slow  convergence,  though,  it  is 
suitable  only  at  singular  points. 

C.  SEMI-INVERSE  INTERACTION  METHOD 

The  semi-inverse  interaction  method  incorporates  direct 
inviscid  and  inverse  viscous  flow  solvers  such  that 
displacement  thickness  is  input  to  both  solvers  as  shown  in 
Figure  6.1c.  External  velocity  is  output  from  both  solvers. 
Convergence  is  ensured  with  the  use  of  a  relaxation  formula 
which  redefines  the  displacement  thickness  distribution. 

Ucv(x) 

6ne„*(x)  =  60i<j*(x)  [1  +  U( -  -  1)  (6.1) 

Uex (x) 

where  u  is  a  relaxation  parameter. 

The  numeric..  1  weaknesses  of  the  direct  and  inverse 
methods  are  improved,  but  inviscid  and  viscous  regions  are 
still  loosely  coupled. 

D.  VISCOUS-INVISCID  INTERACTION  METHOD 

The  viscous-inviscid  interaction  method  ensures  a  strong 
interaction  between  the  outer,  inviscid,  and  inner,  viscous, 
regions.  Both  the  external  velocity,  ue(x),  and  displacement 
thickness,  6*,  are  unknown  quantities.  Convergence  is 
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ensured  through  the. interaction  law  which  uses  the  blowing 
velocity  concept. 

The  equations  are  solved  through  successive  sweeps  over 
the  airfoil  surface  as  indicated  in  Figure  6.2.  For  each 
sweep  the  external  velocity  for  the  boundary  layer  equation 
is  written 


u«(x)  =  u^itx)  +  du^x) 


(6.2) 


where , 


u«x(x)  is  the  inviscid  velocity 


and 


u«(x)  is  the  perturbation  velocity  due  to  the 
boundary  layer  displacement. 

The  perturbation  velocity  is  modeled  by  the  interaction 
law  with  the  help  of  blowing  velocities.  The  displacement 
effect  of  ~  boundary  layer  is  obtained  by  ejecting  fluid  at 
the  surface  of  the  airfoil  as  shown  in  Figure  6.3. 

With  a  properly  arranged  blowing  velocity  source 
distribution  on  the  airfoil  surface,  the  virtual  displacement 
body  becomes  a  streamline. 

In  determining  the  source  strengths,  the  displacement- 
body  tangential  flow  condition  is  represented  by 

v(x,6*)  d6* 
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u.l 


u..  <5* 


VISCOUS  FLOW  SOLVER 

Boundary  laysr  aquation*  ara  tolvad  by  a  finite  difference  method. 


DIRECT  B.L  METHOD 

(U*ad  cloia  to  lead,  edgal 

INPUT:  External  velocity, 
initial  guest  of 
velocity  profile 

OUTPUT:  Velocity  profile  and 
b.l.  characteriatica 
like  dieplacement 
thickness  <5*  and 
akin  friction 


INTERACTIVE  B.L  METHOD 

(Used  anywhere  else) 

Bath  external  velocity  and  displacement 
thickness  are  treated  as  unknowns. 

INPUT:  External  velocity  distribution 
end  disolecement  thickness 
distribution  of  current  (uostreaml 
and  previous  idownstreemi  cycle, 
initial  guess  of  velocity  protile 

OUTPUT:  Velocity  profile  including  tha 
'viscous'  sxtsrnsl  velocity  u« 
and  b.l.  characteristics  like  dispf. 
thickness  4'  and  skin  friction 


OUTPUT  of 

1.  Pressure  distribution 
2.  B.l.  characteristics  for  both 
upper  and  lowar  surfscs 
3.  Velocity  protilea  if  required 


C  STOP  ) 


Figure  6.2  Viscid/Inviscid  Interaction  Method 
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Model  simplifications  are  achieved  through  the  use  of 
these  thin  airfoil  approximations: 

1)  The  u-velocity  component  is  considered  invariant  across 
the  boundary  layer,  as  the  displacement  thickness  is 
thin  enough  to  consider  the  differences  negligible. 

2)  The  blowing  velocity,  v(x,0),  is  one  half  the  source 
strength,  a{x),  on  an  airfoil  represented  by  a  straight 
line . 


o(x) 

-  =  ( x ,  0 ) 

2 


=  v ( x , 6* ) 


‘6*  6v 
—  dy 
o  6y 


d6*  du« 

=  Ue -  +  -  6* 

6x  dx 


a(x)  d 

-  =  - {  Ue6*  ) 

2  dx 


(6.4) 


where  (d/dx)(u«6*)  is  the  blowing  velocity. 

The  blowing  velocity  once  obtained  from  the  source 
strength  is  then  related  to  the  perturbation  velocity,  6u„, 
through  the  use  of  the  Hilbert  integral. 


6u„  =  — 
2rt  J 


°U) 

x-5 


(6.5) 


After  substituting  equations  (6.4)  and  (6.5)  into  (6.2) 
the  interaction  law  is  obtained. 
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(6.6) 


1  fxb  d  d£ 

u«(x)  =  Uei(x)  +  -  — ( ue6* ) - 

n  Jx*  d^  x-J, 

The  numerical  implementation  of  the  interaction  law 
requires  some  discrete  approximation  of  the  thin  airfoil 
integral,  equation  (6.6).  Similar  to  the  panel  method,  a 
piecewise  approximation  of  the  continuous  blowing  velocity 
d(u„6*)/dx  allows  for  piecewise  analytical  integration. 


63 


VII.  AIRFOIL  STUDIES 


Cebeci's  interactive  aircode  was  applied  to  four  airfoils 
over  a  wide  range  of  Reynolds  numbers  and  angles  of  attack. 
The  computer  program  results  were  then  compared  to  reported 
experimental  data.  Unless  otherwise  stated,  20  iterations 
were  used  for  each  computer  run,  and  laminar-to-turbulent 
transitions  were  determined  internal  to  the  program.  The 
significance  of  the  number  of  iterations  will  be  discussed 
later  in  the  section. 

A.  NACA  663-018 

Computer  results  of  the  NACA  663-018  airfoil  section  were 
compared  to  the  test  results  of  Gault  [Ref.  8],  which  were 
performed  in  the  NASA  Ames  Research  Center  7-by-10-foot  wind 
tunnel.  The  laminated  pine  model  with  a  1/8  inch-thick 
mahogany  plywood  veneer  spanned  the  7-foot  dimension  to 
simulate  two-dimensional  flow. 

Total-and  static-pressure  surveys,  hot-wire-anemometer 
observations,  and  detailed  pressure-distribution  and  liquid- 
film  measurements  were  made  in  regions  of  separated  flow. 
The  measurements  were  obtained  for  a  wide  range  of  angles  of 
attack  and  for  Reynolds  numbers  from  1.5  to  10  million.  A 
main  purpose  of  these  measurements  was  to  identify  locations 
of  separation,  transition  and  reattachment. 


64 


Using  the  Cebeci  Code  the  663-018  airfoil  shown  in  Figure 
7.1  was  initially  tested  for  section  lift  coefficients. 
Comparisons  were  made  with  Abbott  and  Doenhoff  [Ref.  9]  and 
the  results  are  presented  in  Figures  7.2  and  7.3  for  Reynolds 
numbers  of  3  and  6  million,  respectively. 

Upper  surface,  laminar  to  turbulent,  transition  locations 
are  shown  in  Figures  7.4  and  7.5  for  increasing  angles  of 
attack  and  for  Reynolds  numbers  of  3  and  6  million, 
respectively.  Gault's  locations  were  obtained  from  pressure 
and  hot-wire  measurements,  which  provided  near  identical 
results.  The  program  transition  locations  were  computed  to 
be  the  point  of  laminar  separation.  Note  that  the  transition 
locations  shift  forward  as  the  angle  of  attack  is  increased, 
and  they  approach  the  leading  edge  above  6  degrees.  Unless 
otherwise  stated,  all  computer  runs  used  a  transition 
constant  of  GY  =  1200. 

Midchord  upper  surface  transitions  at  less  than  two 
degrees  angle  of  attack  and  Reynolds  numbers  of  1.5  to  10 
million  are  shown  in  Figures  7.6,  7.7,  7.8  and  7.9.  In  all 
cases  the  computed  predictions  were  forward  of  Gault's 
because  of  laminar  separation  predicted  by  the  Code. 

while  Gault  found  leading  edge  separation  bubbles,  the 
Cebeci  Code  did  not  predict  them  at  any  angle  of  attack  for 
Reynolds  numbers  of  3  and  6  million. 

The  relationships  between  separation  and  transition  are 
illustrated  in  Figure  7.10  for  the  results  of  Gault  and  the 
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LIFT  COEFFICIENT 


Lift  Coefficient,  NACA  663-018,  R  =  3  Million 
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Figure  7.3  Lift  Coefficient,  NACA  663-018,  R  =  6  Million 
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NACA  663-018,  R  =  3.0  MILLION 


Figure  7.4  Upper  Surface  Laminar  to  Turbulent 

Transition,  NACA  663-018,  R  =  3  Million 


Figure  7.5  Upper  Surface  Laminar  to  Turbulent 

Transition,  NACA  563~013,  R  =  6  Million 
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Figure  7.6  Upper  Surface  Transition,  Midchord,  NACA 
663-OI8,  R  =  1.5  and  2  Million 
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NACA  663-018  MIDCHORD 


Figure  7.7  Upper  Surface  Transition,  Midchord,  NACA 
663-OI8,  R  =  3  and  4  Million 
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UPPER  SURFACE  TRANSITION 


x/c 

NACA  663-018  MIDCHORD 

Figure  7.8  Upper  Surface  Transition,  Midchord,  NACA 
663-JI8,  R  =  6  and  8  Million 
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Figure  7.9  Upper  Surface  Transition,  Midchord,  NACA 
663-OI8,  R  =  10  Million 
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Figure  7.10  Upper  Surface  Transition  and  Separation, 
Midchord,  NACA  663-013,  R  =  2  Million 
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Cebeci  code.  The  experimental  results  show  separation  prior 
to  transition,  whereas  the  computer  results  predict 
separation  after  transition.  The  importance  of  this 
difference  manifests  itself  in  the  difference  between  the 
measured  and  computed  midchord  bubbles,  pressure 
distributions  and  velocity  profiles  shown  in  Figures  7.11  to 
7.31. 

Midchord  separation  bubbles  for  angles  of  attack  of  0  and 
2  degrees  and  Reynolds  number  of  2  million  are  plotted  in 
Figures  7.11  and  7.12.  The  lines  represent  contours  where 
u/u«  =  0.  The  Cebeci  Code  midchord  bubbles  are  much  smaller 
than  those  found  by  Gault. 

Full  chord  pressure  distributions  for  angles  of  attack  of 
zero  and  two  degrees,  and  Reynolds  numbers  of  three  and  six 
million  are  shown  in  Figures  7.13  to  7.16.  In  each  case  the 
biggest  difference  between  the  experimental  results  and  the 
Cebeci  code  occurred  near  the  midchord  separation  bubble 
regions.  Figure  7.17  shows  a  leading  edge  pressure 
distribution  for  an  angle  of  attack  of  six  degrees  and  a 
Reynolds  number  of  three  million. 

Midchord  velocity  profiles  are  shown  in  Figures  7.18 
through  7.24  for  an  angle  of  attack  of  zero  degrees  and  a 
Reynolds  number  of  two  million,  and  in  Figures  7.25  through 
7.31  for  two  degrees  angle  of  attack  and  a  Reynolds  number  of 
two  million.  These  velocity  profiles  clearly  show  a  big 
difference  in  bubble  sizing. 
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MIDCHORD  BUBBLE  SHAPE 
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NACA  663-018,  A0A=0,  R=2  MILLION 


Figure  7.11  Midchord  Bubble  Shape,  NACA  663 -018, 
AOA  =  0°,  R  =  2  Million 
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NACA  663-018,  A0A=2,  R=2  MILLION 


Figure  7.12  Midchord  Bubble  Shape,  NACA  663-OI8, 
AOA  =2°,  R  =  2  Million 
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x/c 

NACA  663-018  A0A=0  R=3  MILLION 


Figure  7.13  Upper  Surface  Pressure  Distribution,  NACA 
663-OI8,  AOA  =  0°,  R  =  3  Million 
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UPPER  SURFACE  PRESSURE 


NACA  663-018  AOA=0  R=6  MILLION 


Figure  7.14  Upper  Surface  Pressure  Distribution,  NACA 
663-018,  AOA  =  0°,  R  =  6  Million 
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UPPER  SURFACE  PRESSURE 


x/c 

NACA  663-018  A0A=2  R=3  MILLION 


Upper  Surface  Pressure  Distribution,  NACA 
66*3-018,  AOA  =2°,  R  =  3  Million 
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NACA  663-018  A0A=2  R=6  MILLION 


Figure  7.16  Upper  Surface  Pressure  Distribution,  NACA 
663-OI8,  AOA  =  2°,  R  =  6  Million 
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NACA  663-018  A0A=6  R=3  MILLION 


Figure  7.17  Leading  Edge  Upper  Surface  Pressure 

Distribution,  NACA  663-018,  AOA  =  6°, 
R  =  3  Million 
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Figure  7.20  Upper  Surface  Velocity  Profil 

X/C  =  .64,  AOA  =0°,  R  =  2  Million 
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VELOCITY  RATIO,  U/UE 

NACA  663-018,  AOA=0,  R=2  MILLION 


Figure  7.22 


Upper  Surface  Velocity  Profile,  NACA  663-018, 
X/C  =  .68,  AOA  =0°,  R  =  2  Million 


VELOCITY  PROFILES 


VELOCITY  RATIO,  U/VE 

NACA  663-018,  A0A=2,  R=2  MILLION 


Upper  Surface  Velocity  Profile,  NACA  66 
X/C  =  .58,  ACA  =  2°,  R  =  2  Million 
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VELOCITY  RATIO,  U/l'E 


NACA  663-018,  A0A=2,  R=2  MILLION 

Figure  7.23  Upper  Surface  Velocity  Profile,  NACA  66 
X/C  =  .64,  AOA  =  2°,  R  =  2  Million 
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VELOCITY  RATIO,  U/UE 

NACA  663-018,  A0A=2,  R=2  MILLION 


Figure  7.29  Upper  Surface  Velocity  Profile,  NACA  663-C18, 
X/C  =  .66,  AOA  =  2°,  R  =  2  Million 
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VELOCITY  RATIO,  UA'E 

NACA  663-018,  AOA-2,  R=2  MILLION 


Figure  7.30  Upper  Surface  Velocity  Profile,  NACA  663 
X/C  =  .68,  AOA  =  2°,  R  =  2  Million 


In  an  effort  to  increase  the  region  of  separated  flow, 
both  the  GY  transition  constant  and  the  location  of  upper 
surface  transition  were  adjusted.  Table  7.1  shows 
theresults.  With  a  constant  of  1200  the  transition  location 
could  not  be  moved  aft.  However,  when  the  constant  was 
lowered  to  200  and  below,  the  transition  location,  x/c  =  .69, 
found  by  Gault,  could  be  used  moving  the  transition  inside 
the  bubble.  While  a  lowered  GYtt-  constant  and  increased  x/c 
transition  improved  the  bubble  size,  the  separation  length, 
x/c  =  .60  to  .70,  could  not  quite  be  met.  The  best  result, 
bubble  =  .6391  -  .7098  (x/c),  was  obtained  with  a  GYtr-  of  40, 
and  an  upper  surface  transition  location,  XTRU,  input  of  .69 
(x/c).  Whether  40  is  a  suitable  value  for  other  foils  has 
not  been  determined. 

Twenty  iterations  were  used  on  all  computer  runs  for  this 
airfoil.  To  make  sure  that  20  iterations  were  sufficient  for 
accurate  results.  Figure  7.32  was  obtained.  The  lift 
coefficient  was  plotted  for  each  iteration,  and  as  can  be 
seen  the  even  iterations  produced  very  minimal  changes  past 
12.  Even  between  19  and  20  the  change  in  lift  was  only  5  x 
10-*.  Therefore,  20  iterations  were  considered  sufficient 
for  all  computer  runs. 

B.  NACA  0010  (MODIFIED) 

Similar  to  the  NACA  66s-018,  computer  results  of  the  NACA 
0010  (Modified)  airfoil  section  were  compared  to  the  test 
results  of  Gault  [Ref.  8J.  The  tests  were  also  conducted  in 
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TABLE  7 . 1  EFFECT  OF  GYTR  AND  XTRU  ON  THE  LENGTH  OF  THE 
SEPARATION  BUBBLE 


GYTR 

XTRU  (X/C) 

SEPARATION 

(X/C) 

1200 

.620792 

(1) 

.6572  -  . 

6929 

.65 

(2) 

.69 

(2) 

400 

.69 

(2) 

300 

.  69 

(2) 

200 

.69 

.6391  - 

.7596 

120 

.639164 

(1) 

No  Separation 

.65 

.6572  - 

.6751 

.69 

.6391  - 

.7434 

80 

.69 

.6391  - 

.7268 

60 

.69 

.6391  - 

.7268 

40 

.69 

.6391  - 

.7098 

30 

.69 

.6572  - 

.7098 

20 

.69 

.6572  - 

.7098 

10 

.69 

.6572  - 

.7098 

NACA  663-018 

Reynolds  Number  =  2,000,000 
Angle  of  Attack  =  0  Degrees 

Region  of  Separation,  Experimental  (Gault)  =  .60  -  .70  (x/c) 
GYTR  =  Empirical  Transition  Constant 

XTRU  =  Upper  Surface  Begin  of  Transition,  Input  or  Computer 
Derived 

Lower  Surface  Begin  of  Transition  =  Computer  Derived,  Each 
Case 

1 )  Computer  Derived 

2)  Breakdown  in  Simulation 
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Figure  7.32  Lift  Coefficient  Versus  Iterations,  NACA 

663-OI8,  AOA  =  0,  R  =  2  Million,  Transition 
Constant  =  1200 
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the  7-by-10-foot  wind  tunnel  at  NASA  Ames  and  two-dimensional 
flow  was  simulated. 

The  0010  (modified)  airfoil  shown  in  Figure  7.33,  unlike 
the  663-018,  was  not  computer  tested  for  section  lift 
coefficients  as  Abbott  and  Doenhoff  [Ref.  9]  had  no 
comparable  airfoil. 

Leading  edge,  upper  surface,  laminar  to  turbulent, 
transition  locations  were  observed,  though,  and  the  results 
are  shown  in  Figures  7.34  and  7.35  for  Reynolds  numbers  of 
two  and  six  million,  respectively.  The  Cebeci  Code  curves 
represent  the  beginning  of  transition. 

Full  chord  pressure  distributions  for  angles  of  attack  of 
zero  and  three  degrees  and  Reynolds  numbers  of  three  and 
eight  million  are  plotted  in  Figures  7.36  through  7.39. 

Leading  edge  pressure  distributions  for  angles  of  attack 
of  four,  eight  and  twelve  degrees,  and  Reynolds  numbers  of 
two  and  six  million  are  shown  in  Figures  7.40  through  7.45. 
Of  particular  interest  are  the  "lump"  disparities  in  Figures 
7.42  through  7.45.  A  possible  explanation  for  the  computer 
program  deletions  of  the  lumps  is  a  failure  to  predict 
leading  edge  bubbles. 

C.  NACA  4412 

Computer  results  of  the  NACA  4412  airfoil  section  were 
compared  to  the  test  results  of  Hastings  and  Williams  [Ref. 
10],  which  were  performed  in  the  13-by  9-foot  low  speed  wind 
tunnel  of  the  Royal  Aircraft  Establishment  at  Bedford.  The 
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LAMINAR  TO  TURBULENT  TRANSITION 


Figure  7.34  Upper  Surface  Laminar  to  Turbulent 
Transition,  Leading  Edge,  NACA  0010 
(Modified) ,  R  =  2  Million 
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NACA  0010  MOD,  R  =  6.0  MILLION 


Figure  7.35  Upper  Surface  Laminar  to  Turbulent 
Transition,  Leading  Edge,  NACA  0010 
(Modified),  R  =  6  Million 
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NACA  0010  (MOD)  A0A=0  R=3  MIL 


Figure  7.36  Upper  Surface  Pressure  Distribution,  NAC 
0010  (Modified),  AOA  =  0°,  R  =  3  Million 


Figure  7.37  Upper  Surface  Pressure  Distribution,  NACA 
0010  (Modified),  AOA  =  0°,  R  =  8  Million 


105 


UPPER  SURFACE  PRESSURE 
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NACA  0010  (MOD)  A0A=3  R=8  MIL 

figure  7.39  Upper  Surface  Pressure  Distribution,  NACA 
0010  (Modified),  AOA  =  3°,  R  =  3  Million 
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PRESSURE  COEFFICIENT 
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UPPER  SURFACE  PRESSURE 
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NACA  0010  (MOD)  A0A=4  R=2  MIL 

Figure  7.40  Leading  Edge  Upper  Surface  Pressure 

Distribution,  NACA  0010  (Modified),  AOA  =  4°, 
R  =  2  Million 
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NACA  0010  (MOD)  A0A=4  R=6  MIL 


Figure  7.41  Leading  Edge  Upper  Surface  Pressure 

Distribution,  NACA  0010  (Modified),  AOA  =  4°, 
R  =  6  Million 
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NACA  0010  (MOD)  AOA=8  R=2  MIL 


Figure  7.42  Leading  Edge  Upper  Surface  Pressure 

Distribution,  NACA  0010  (Modified),  AOA 
R  =  2  Million 


NACA  0010  (MOD)  A0A=8  R=6  MIL 


Figure  7.43  Leading  Edge  Upper  Surface  Pressure 

Distribution,  NACA  0010  (Modified),  AOA 
R  =  6  Million 


UPPER  SURFACE  PRESSURE 


Figure  7.44  Leading  Edge  Upper  Surface  Pressure 
Distribution,  NACA  0010  (Modified), 
AOA  =  12°,  R  =  2  Million 
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Figure  7.45  Leading  Edge  Upper  Surface  Pressure 
Distribution,  NACA  0010  (Modified), 
AOA  =  12°,  R  =  6  Million 


113 


one  meter  chord  model  spanned  the  full  width,  13  feet,  to 
simulate  two-dimensional  flow. 

Mounted  at  its  quarter-chord  point,  the  model  was 
extensively  instrumented  with  static  pressure  orifices. 
Boundary  layer  and  wake  measurements  were  made  at  mid-span 
where  the  88  pressure  orifices  were  located. 

The  main  emphasis  in  the  experiment  was  on  defining  the 
upper  surface  boundary  layer  through  separation  and  into  the 
wake.  Laser  anemometry  was  used  to  measure  the  average 
velocities  and  Reynolds  stresses  were  measured  by  hot-wire 
anemometry. 

The  4412  shown  in  Figure  7.46  was  initially  computer 
tested  with  the  Cebeci  code  for  momentum  thickness.  In 
Figure  7.47  the  upper  and  lower  surface  laminar  to  turbulent 
transitions  were  computer  derived,  x/ctjr  upper  and  lower 
surfaces  =  0.00625.  In  Figure  7.48  the  upper  and  lower 
surface  transitions  were  input  as  x/ctr-  upper  surface  =  0.01, 
and  x/Ctr  lower  surface  =  0.11.  These  values  were  as  close 
as  could  be  input  to  0.014  and  0.110  respectively,  for  the 
downstream  ends  of  the  transition  trips  used  in  the 
experiment.  The  differences  in  transition  locations  seemed 
to  make  no  difference  in  computer  results.  The  momentum 
thicknesses  still  did  not  agree  very  well  with  Hastings  and 
Williams'  experimental  results. 

Figure  7.49  compares  lift  coefficients  from  the  Cebeci 
code,  Hastings  and  Williams,  and  Abbott  and  Doenhoff  [Ref. 
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Figure  7.47  Upper  Surface  Momentum  Thickness,  NACA  4412, 
AOA  =  12.49°,  R  =  4.17  Million,  Computer 
Derived  Transitions 


115 


THKTA/C  X1000 

2  0  3.0  4.0  SO 


UPPER  SURFACE  MOMENTUM  THICKNESS 


x/c 

NACA  4412  R=4.17  MIL  AOA=  12.49 

Figure  7.48  Upper  Surface  Momentum  Thickness,  NACA  4412, 
AOA  =  12.49°,  R  =  4.17  Million,  Computer 
Derived  Input  Transitions 
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Figure  7.49  Lift  Coefficients.  NACA  4412 


10].  The  dashed  line  for  the  Cebeci  code,  Reynolds  number 
equal  to  4.17  million  and  computer  derived  transitions, 
should  lie  between  the  curves  for  Reynolds  number  equal  to 
three  and  six  million  from  Abbott  and  Doenhoff.  However,  it 
lies  above  the  six  million  curve  which  further  indicates  an 
insufficient  boundary  layer  development.  With  the  input 
transition  the  code  prediction  for  Reynolds  number  equal  to 
4.17  million  and  angle  of  attack  equal  to  12.49  degrees,  was 
even  higher.  Of  interesting  note  though,  is  that  the 
Hastings  and  Williams  prediction,  Reynolds  number  equals  4.17 
million  and  angle  of  attack  equals  12.49  degrees,  lies  below 
the  Abbott  and  Doenhoff  values,  possibly  indicating  an  error 
on  their  part. 

Figures  7.50  through  7.56  show  the  velocity  profiles  from 
x/c  =  .66  to  the  trailing  edge.  In  all  cases  the  Reynolds 
number  was  4.17  million,  the  angle  of  attack  was  12.49 
degrees,  and  the  upper  and  lower  transitions  were  .01  and 
.11,  respectively.  U/Ue  indicates  the  fraction  of  the 
velocity  at  the  boundary  layer  edge,  and  ri /Delta  is  the 
fraction  of  boundary  layer  thickness,  where  Delta,  6,  is 
defined  as  the  layer  thickness  where  the  velocity  is  99%  of 
the  edge  velocity. 

As  these  figures  indicate,  as  well  as  Figure  7.57,  a 
Cebeci  code  velocity  profile  summation,  the  code  does  predict 
separation,  but  not  the  extent  indicated  by  Hastings  and 
Williams.  If  the  lift  coefficient  curves,  Figure  7.58,  can 
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ETA/DELTA 

NACA  4412  R=4.17  MIL  AOA=12.49 

Figure  7.50  Upper  Surface  Velocity  Profile,  NACA  4412, 
X/C  =  .66,  AOA  =  12.49°,  R  =  4.17  Million 
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0  0 


ETA/DELTA 

NACA  4412  R=4.17  MIL  A0A=12.49 


Figure  7.52  Upper  Surface  Velocity  Profile,  NACA  4412, 
X/C  =  .78,  ACA  =  12.49°,  R  =  4.17  Million 
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UPPER  SURFACE  VELOCITY  PROFILE 


NACA  4412  R=4.17  MIL  A0A=12.49 


Figure  7.53  Upper  Surface  Velocity  Profile,  NACA  4413, 
X/C  =  .85,  AOA  =  12.49°,  R  =  4.17  Million 
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Figure  7.54  Upper  Surface  Velocity  Profile,  NACA  4412, 
X/C  =  .92,  AOA  =  12.49°,  R  =  4.17  Million 
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NACA  4412  R=4.17  MIL  AOA=12.49 


Figure  7.55  Upper  Surface  Velocity  Profile,  NACA  4412, 
X/C  =  .95,  AOA  =  12.49°,  R  =  4.17  Million 
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UPPER  SURFACE  VELOCITY  PROFILE 


Figure  7.56  Upper  Surface  Velocity  Profile,  NACA  4412, 
X/C  =  .997,  AOA  =  12.49°,  R  =  4.17  Million 
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NACA  4412,  CEBECI  CODE 


Figure  7.57  Upper  Surface  Velocity  Profiles,  NACA  4412, 
AOA  =  12.49°,  R  =  4.17  Million 
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LIFT  COEFFICIENT 


co 


ANGLE  OF  ATTACK 

NACA  4412  R  =  1.523  MILLION 


Figure  7.58  Lift  Coefficient,  NACA  4412,  20  Iterations, 
R  =  1.523  Million 
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be  a  reference,  then  it  appears  that  the  Cebeci  code  predicts 
an  underdeveloped  flow,  too  little  separation,  and  Hastings 
and  Williams  show  an  overdeveloped  flow,  too  much  separation, 
for  the  given  conditions. 

To  insure  that  the  Cebeci  Code  was  run  correctly,  certain 
results  by  Cebeci,  Clark,  Chang,  Halsey  and  Lee  [Ref.  1]  were 
attempted  to  be  duplicated.  Figure  7.58  compares  two  curves 
from  Figure  14,  [Ref.  1],  curves  labeled  interactive  theory 
and  interactive  theory  with  a  modified  transition,  with  a 
curve  obtained  using  the  Cebeci  Code  with  20  iterations  and 
computer  derived  transitions.  Interestingly,  the  first  and 
third  curves  of  Figure  7.58,  interactive  theory  and  Cebeci 
Code,  respectively,  should  be  the  same,  but  the  two  clearly 
are  not  above  nine  degrees  angle  of  attack.  Even  more 
interestingly,  the  Cebeci  Code  lift  curve  in  Figure  7.59 
after  only  10  iterations  does  match  the  interactive  theory 
curve . 

Figure  7.60  clearly  shows  the  importance  of  using  enough 
iterations  to  obtain  a  reasonably  accurate  solution. 

Finally,  Figure  7.61  shows  a  very  good  match  between  the 
pressure  coefficients  for  set  conditions  of  Figure  16,  Cebeci 
et  al  [Ref.  1],  and  the  Cebeci  Code,  20  iterations. 

D.  FX  63-137 

Computer  results  of  the  Wortmann  FX  63-137  airfoil  were 
compared  to  the  test  results  of  Brendel  and  Mueller  [Ref. 
11],  which  were  conducted  in  the  University  of  Notre  Dame 
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LIFT  COEFFICIENT 
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Figure  7.59  Lift  Coefficient,  NACA  4412,  20  Iterations, 
R  =  1.523  Million 
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ITERATIONS 

NACA  4412  R=1.523  MILLION  A0A=12 


Figure  7.60  Lift  Coefficient  Versus  Iterations, 

NACA  4412,  AOA  =  12°,  R  =  1.523  Million 
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.61m  x  .61m  wind  tunnel.  Two  cast  epoxy  resin  airfoil  models 
with  chords  of  .305m  and  spans  of  .4m  were  mounted  in  the 
center  of  the  test  section.  Pressure  was  recorded  on  one 
model  with  96  pressure  taps  connected  through  two  scanivalves 
to  an  electronic  manometer.  Boundary  layer  velocity 
measurements  were  obtained  on  the  other  model  using  a 
constant  temperature  anemometer  with  a  five  pm  diameter, 
single-sensor,  hot-wire,  boundary  layer  probe. 

Using  the  Cebeci  Code  the  FX  63-137  airfoil  shown  in 
Figure  7.62  was  initially  tested  for  section  lift 
coefficients  with  a  transition  constant  of  1200.  Reynolds 
numbers  of  .28,  .5  and  .7  million  were  used,  and  the  results, 
shown  in  Figures  7.63  and  7.64,  were  compared  to  those  of 
Aithaus  and  Wortmann  [Ref.  12].  Interestingly,  the  Cebeci 
Code  predicted  low  values  for  Reynolds  numbers  of  .28  and  .5 
million,  but  for  .7  million  the  lift  coefficients  were  nearly 
identical  to  Aithaus  and  Wortmann  up  to  an  angle  of  attack  of 
10  degrees. 

As  the  purpose  of  Brendel  and  Mueller  was  to  make 
boundary  layer  measurements  at  low  Reynolds  numbers,  a 
computer  comparison  was  unsuccessfully  attempted  for  steady 
flow  at  a  Reynolds  number  of  100,000  and  an  angle  of  attack 
of  7  degrees.  With  20  iterations  the  Cebeci  Code  failed. 

To  understand  why  the  code  calculations  ceased  for  this 
case,  other  computer  runs  were  attempted  for  the  same 
Reynolds  number  and  angle  of  attack,  but  with  fewer 
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Figure  7.63  Lift  Coefficient,  FX  63-137,  R  =  .28  Million, 
and  .7  Million 
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Figure  7.64  Lift  Coefficient,  FX  63-137,  R  =  .5  Million 
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iterations.  Figures  7.65  and  7.66  show  the  upper  surface 
displacement  and  momentum  thicknesses  for  steady  flow  and 
iterations  of  2,4,  6,  8  and  10.  As  can  be  seen  in  both 
figures  flow  calculations  matched  very  well  with  experimental 
data  up  to  approximately  x/C  =  .55.  After  that  point  stall 
occurred  and  calculations  ceased  with  more  than  10 
iterations.  Brendel  and  Mueller  experimentally  derived 
separation  to  begin  at  x/C  =  .34,  but  reattachment  was  shown 
to  occur  at  x/c  =  .60.  Unfortunately  the  Cebeci  Code  could 
not  predict  reattachment  for  the  prescribed  conditions. 
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Figure  7.65  Upper  Surface  Displacement  Thickness, 
FX  63-137,  AOA  =  7°,  R  =  100,000 
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UPPER  SURFACE  MOMENTUM  THICKNESS 


;  ; 1 l  i - ; :  r  :  :  ; 

0.0  0.1  0.2  0.3  0.4  0.5  0.8  0.7  0.8  0.9  1.0 

X/C 

FX  63-137  R  =  100,000  AOA  =  7 


Figure  7.66  Upper  Surface  Momentum  Thickness,  FX  63-137, 
AOA  =  7° ,  R  =  100,000 
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VIII.  CONCLUSIONS  AND  RECOMMENDATIONS 


Cebeci's  viscous/inviscid  interaction  program  was  applied 
to  the  analysis  of  steady,  two  dimensional,  incompressible 
flow  past  four  airfoils,  the  NACA  663-018,  0010  (modified), 
4412  and  the  Wortmann  FX  63-137.  Detailed  comparisons  with 
the  available  experimental  results  show  that  for  attached 
flows  the  essential  features  are  correctly  modelled,  but  that 
significant  discrepancies  are  found  in  regions  of  flow 
separation.  These  discrepancies  are  possibly  caused  by  the 
empirical  transition  modelling  used  in  the  present  code. 
Future  efforts  therefore  should  be  directed  to  the 
incorporation  of  transition  calculations  which  permit  the 
prediction  of  transition  within  a  separation  bubble,  such  as 
the  application  of  the  e"-method  proposed  by  Cebeci  [Ref. 
13] . 
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